Hybridization may aid evolutionary rescue of an endangered East African passerine

Abstract Introgressive hybridization is a process that enables gene flow across species barriers through the backcrossing of hybrids into a parent population. This may make genetic material, potentially including relevant environmental adaptations, rapidly available in a gene pool. Consequently, it has been postulated to be an important mechanism for enabling evolutionary rescue, that is the recovery of threatened populations through rapid evolutionary adaptation to novel environments. However, predicting the likelihood of such evolutionary rescue for individual species remains challenging. Here, we use the example of Zosterops silvanus, an endangered East African highland bird species suffering from severe habitat loss and fragmentation, to investigate whether hybridization with its congener Zosterops flavilateralis might enable evolutionary rescue of its Taita Hills population. To do so, we employ an empirically parameterized individual‐based model to simulate the species' behaviour, physiology and genetics. We test the population's response to different assumptions of mating behaviour and multiple scenarios of habitat change. We show that as long as hybridization does take place, evolutionary rescue of Z. silvanus is likely. Intermediate hybridization rates enable the greatest long‐term population growth, due to trade‐offs between adaptive and maladaptive introgressed alleles. Habitat change did not have a strong effect on population growth rates, as Z. silvanus is a strong disperser and landscape configuration is therefore not the limiting factor for hybridization. Our results show that targeted gene flow may be a promising avenue to help accelerate the adaptation of endangered species to novel environments, and demonstrate how to combine empirical research and mechanistic modelling to deliver species‐specific predictions for conservation planning.

Individuals refer to individual birds belonging to either of the species Z. silvanus and Z. avilateralis. Each individual has a genome consisting of multiple genes, which each code for a specic trait. Genes are grouped into linkage units (chromosomes), although for the purposes of most experiments in this study we congured the model to create a separate linkage unit for each gene. Individuals are diploid, so their phenotype for each trait is determined by the mean value of the two genes coding for this trait. These traits determine an individuals adaptation to its habitat (i.e. temperature and above-ground carbon) and inuence life-history processes like dispersal and reproduction (Fig. S1). Individuals are assigned a sex at birth and may nd a partner after dispersal.
The model simulates the landscape of the Taita Hills (total area: 962 km², or 96,170 patches) over 300 years, with one time step per year.  Figure S1: Relationship of genotype and phenotype in GeMM individuals. Note that the model was used in a simplied setup in this study, so not all of the mentioned traits are relevantsee text for details. (Modied from Leidinger et al., 2021.) 1.3 Process overview and scheduling Patches and their communities are created during initialisation. The patch environment stays constant, but each year, each patch community undergoes the following processes (in order): 1. Establishment: Individuals that dispersed in the previous year are checked for genetic viability. Any individuals that are not viable due to extreme low tness or genetic defects are killed.
2. Survival: Density-independent mortality of individuals. 3. Reproduction: Each breeding pair in a patch produces up to two ospring. Parent genomes undergo meiosis, leading to recombination in the ospring genomes. 5. Dispersal: Juveniles leave the parent patch, iteratively moving to the neighbouring patch whose habitat is closest to their optimum. If they nd a patch with an available mate or a free territory, they join the patch community. If they reach their maximum dispersal distance without settling down, they die.

Design concepts
Basic principles.
To estimate unknown parameters and process rates, GeMM uses the Metabolic Theory of Ecology (MTE; Brown et al., 2004). This describes biological rates such as reproduction, mortality, and mutation as a mathematical function of temperature and individual body mass. As we have good data for most Zosterops life processes, and to avoid additional complexity, we only applied the MTE to the mutation submodel in this study. Emergence. The spatial and temporal demographics and population genetics of the two Zosterops species are emergent properties of the model, arising from the interaction of individuals with each other and the landscape. Objectives.
During dispersal, individuals look for patches with suitable habitats based on their above-ground carbon (AGC) preference. Goodness of t is dened as the difference between patch AGC and individual optimum AGC, with smaller values being better. Sensing.
During dispersal, individuals perceive the AGC values of each neighbouring patch (i.e. 100 m sight radius). Within the current patch, individuals sense the availability of free space and the suitability of another individual as a partner.

Interactions.
A male and a female individual can form a breeding pair, exclusively mating with each other while both are alive. Indirect interaction takes place in the form of competition for space and partners, as individuals that fail to nd an available territory or partner die.

Stochasticity.
Random components inuence the initial population sizes, the likelihood of per-individual density-independent mortality, the reduction of the diploid genomes into haploid gametes during meiosis, and the number of ospring per breeding pair per year. Observation.
The model records per-patch population sizes over time, as well as the distribution of trait values and the degree of heterozygosity in each patch population.

Initialisation
The rst step of model initialisation is to read in all model parameters (see below). These include species denitions for all Zosterops species that are to appear in the model, which are used to create species archetype objects.
Following this comes the process of reading in the user-dened map le. This species the properties of each individual patch, including which patches are to be initialised with a community, and the local carrying capacity. Once the patch object has been created, if it is to have an initial community, the model checks whether its habitat is suitable for any of the predened Zosterops species.
If that is the case, the initial community size for the patch is determined at random. The community is then created by copying the suitable archetypes, always generating an equal number of males and females to form the rst breeding pairs. After copying, the new individuals are mutated to give rise to the initial population variability.

Input data
GeMM takes two external input les for each run: the conguration le with relevant model parameters, and the map le specifying the simulation world. The conguration le can include all parameters listed in Appendix 1.8. Parameters that are not explicitly congured revert to the default, as dened in the source code. Notable parameters include the fertility rate, mean dispersal distances, and hybridisation propensity, as well as the species denitions mentioned above. Also, a seed for the random number generator can be set, allowing exact replicates to be run. (Simulations with identical seed and input and will be identical in output.) The default map le for the Taita Hills (cf. Fig. S2) was generated from remote sensing AGC data as published in Adhikari et al. (2017) and Pellikka et al. (2018). The original data was obtained in 2015 through a combination of airborne laser scanning, satellite imagery, and eld work at a resolution of 32 × 32 m. As we assumed one-hectare patches for GeMM, we upscaled this to 100 × 100 m using cubic convolution resampling in QGis (QGIS.org, 2020). The carrying capacity K (in birds/ha) for each patch was calculated from the patch's above-ground carbon value P AGC as: This equation was chosen to relate the AGC ranges of dierent habitat types (Pellikka et al., 2018) to observed breeding densities of Zosterops in these habitats (J. Engler, personal communication; Mulwa et al., 2007). It gives a density of one breeding pair per hectare in woodland, 23 pairs/ha in exotic forest, and 48 pairs/ha in montane forest.
To account for very low population densities in the savannah (observed as one breeding pair per square kilometer), patches with P AGC < 10 had a 1 % chance of being assigned a carrying capacity of 2. The scenario maps for the habitat experiments were based on the same data as the default map, but edited to reect the scenarios described in Section 2.2. (Table S1). Note that the landscape was not dynamic, i.e. scenario maps stayed constant and there was no habitat change throughout the simulation. (While GeMM is capable of using multiple maps over a single simulation run, but we decided not to use this feature due to the complexities of modelling land-use changesee main text.)

Submodels
All relevant model parameters are listed in Table S2, and will be referred to in the following using block print. For a comprehensive list of parameters (including those not relevant to the study presented here), see Listing 2.  Establishment is carried out for all individuals that are new to a patch following dispersal by testing their genetic viability. This includes ensuring that its genome codes for all necessary traits, and that all trait and adaptation values are within acceptable bounds (e.g. >0). If the individual fails any of these tests, it is removed from the patch community and dies. Survival.
Survival represents density-independent mortality of all adult individuals (i.e. post-dispersal, post-establishment). Each year, every adult has a death probability determined by the global mortality parameter. The value for this was chosen to reect known life expectancies of related Zosterops species (Bird et al., 2020).

Reproduction.
Every adult, established individual with a partner reproduces every year. The number of ospring per breeding pair per year is a randomly chosen integer between 0 and fertility. This value takes into account observed Zosterops clutch sizes and annual breeding attempts as well as nestling and edgling mortality rates. It thus represents the number of juveniles who survive until post-natal dispersal.
For each ospring individual, the parent genomes undergo meiosis, producing two randomly assembled haploid chromosome sets that are combined to form the ospring genome. Ospring are assigned a random sex and the species of their parents. In the case of hybrid ospring (with parents of dierent species), the species of the parent they are phenotypically closer to is chosen. Mutation.
If the mutate parameter is true, each juvenile's genome is mutated after birth. Mutations randomly change individuals' phenotypes and may thus be detrimental, neutral, or benecial, depending on the environment. The number of mutations per individual U (I) is calculated using the Poisson distribution as: where µ is the global parameter mutationrate, which is followed by the metabolic coecient with the activation energy E act , the Boltzman constant k B , and the patch temperature P temp (cf. Appendix 1.8; Brown et al., 2004).
Each mutation is then carried out by choosing a random gene and changing one of its nucleotides. Simultaneously, the trait that it codes for is changed using the normal distribution: where T m is the mutated trait value, T o is the original trait value, and ϕ is the global parameter phylconstr.
Note that mutations were turned o for this study, except during one exploratory scenario (see Section 2.1). Dispersal.
The model simulates post-natal dispersal of all juveniles in the year they are born. Secondary (i.e. adult) dispersal is not simulated.
The dispersal algorithm is a simplied version of the Stochastic Movement Simulator (SMS) developed by Palmer et al. (2011). This algorithm has been found to produce realistic dispersal patterns when tested against real movement data of two other forest bird species in the Taita Hills (Cabanis' greenbul Phyllastrephus cabanisi and whitestarred robin Pogonocichla stellata; Aben et al., 2014).
The maximum dispersal distance of each individual is drawn from a logistic distribution: where I dispmean is the individuals mean dispersal distance and I dispshape is the scale parameter. (For the archetype individuals, these values are initially set to dispmean and dispshape, respectively.) The individual keeps moving until it either nds a suitable patch to settle in, or the maximum dispersal distance is reached, in which case it dies. During each movement step, it evaluates the dierence between neighbouring patches' AGC value and its own AGC optimum, moving to the patch with the least dierence. There is no penalty for moving to patches outside the individual's AGC tolerance range. Patches that have already been visited are ignored.
If the current patch is within its AGC tolerance range, the individual checks for an empty territory (i.e. two empty slots in the patch capacity) or an available mate. An individual already living in that patch is an available mate if it is of the opposite sex, does not have a partner, and there is a free slot in the patch capacity. Individuals rst search for available mates among conspecic individuals. If there are no conspecic available mates and the speciation parameter is set to o, the individual then searches for a non-conspecic mate, accepting a potential mate with a probability determined by the tolerance parameter (i.e. the hybridisation propensity). If a suitable patch is found, the individual is added to the patch community. If there is an available mate, the two individuals are assigned each other as partner, thus forming a new breeding pair. Breeding pairs stay faithful for life, but individuals whose partners die can form new breeding pairs with younger individuals that disperse into their patch.

Implementation details
This appendix provides further details on model parameters and variables.
The full model source code along with the scripts that were used to set up and analyse the experiment may be found on the accompanying CD ROM, or online on Github at https://github.com/CCTB-Ecomods/gemm. The raw data has been archived on the CCTB servers under /storage/ecomod/zosterops/hybridstudy.
Listing 1 shows the full list of state variables for individuals and patches. Listing 2 contains all model parameters and their generic value (for values used in the current study, see Table S2).
Second, although this study deals with evolutionary rescue, we disabled mutations by default to prevent conating eects. Nonetheless, mutations are of course important and ought to be considered. Therefore, to compare the eect size of introgressive hybridisation with that of mutation-based natural selection, we simulated three scenarios in which mutation was turned on. For these, we set the hybridisation propensity to 0, while varying the mutation rate parameter between 0, 3.6 × 10 9 , 3.6 × 10 10 , and 3.6 × 10 11 (giving a maximum of 0, 3, 6, and 21 mutations per individual, respectively; see section 1.7). Mutations in GeMM may be detrimental, neutral, or benecial depending on the individual's environment, and are therefore subject to natural selection.
Third, in the main study we also assumed a simplied genetic architecture with no linkage, i.e. one gene per chromosome (Section 1.2). Again, this was done to prevent conating eects. However, it has been shown that genetic architecture can have a strong eect on a species' adaptive potential (Schiers et al., 2014;Uecker et al., 2015); and a previous study with GeMM showed a clear selective benet of low linkage, due to the increased recombination possibilities of unlinked genes (Leidinger et al., 2021). Hence, to verify that the assumption of no linkage did not qualitatively aect the results of the study, we ran additional scenarios with full linkage (where a single chromosome held all genes of a haploid genome) and random linkage (varying numbers of chromosomes). For these scenarios, mutations were turned o and the hybridisation propensity set to 0.01.
Finally, we twice repeated the habitat experiment as described in the main text, but set the hybridisation propensity to 10 % and 0 %, respectively, rather than 1 %.
Each of these exploratory experiments was run with 10 replicates.

Results
When simulating the hybridisation experiment over 1000 years, population levels plateaued after 300 years, with large overlaps between scenarios with hybridisation. Compared to the 300 year experiment, there were changes in the ranking of scenarios' population sizes, with intermediate hybridisation propensities doing slightly better than high hybridisation propensities. However, the dierences remained small compared to the variability (Fig. S3a, S4). Population heterozygosity, AGC optimum, and AGC tolerance were largely the same as after 300 years, and kept the ranking by hybridisation propensity (Fig. S3b,c,d).
In the linkage experiment, scenarios with linkage (either full or random) were quite distinct from the no linkage scenario. Population growth occurred in all scenarios, but was strongest without linkage. Mean heterozygosity was higher with linkage, as was the AGC optimum. After a small initial rise, the AGC tolerance dropped again in scenarios with linkage, while maintaining its higher-than-initial level without linkage (Fig. S5).
In the mutation experiment, populations of Z. silvanus also grew over time, with more growth at higher mutation rates. Growth was slower than in the hybridisation experiment, but, at the highest mutation rate, reached much higher population levels (Fig.  S6a). As in the hybridisation experiment, AGC optimum dropped steadily; however, AGC tolerance rose in all scenarios (Fig. S6c,d).
The repeat of the habitat experiment with a hybridisation propensity of 10 % yielded qualitatively the same results as with a hybridisation propensity of 1 %, and closely followed the trajectories of the hybridisation experiment at 10 % hybridisation propensity (Fig. S7). At 0 % hybridisation propensity, no intermixing of species took place, and population development consequently mirrored that seen in the hybridisation experiments in the equivalent scenario, although the nal population size reached depended on the amount of montane forest habitat available (Fig. S8).   To allow comparisons from our simulation results to genomic inferences, we performed Pairwise sequential Markovian coalescent (PSMC) analyses (Patton et al., 2019) of historical eective population size (N e ) to test for signatures of potential hybridization between the highland Z. silvanus and its lowland congener Z. avilateralis. We analysed ve individuals of four closely related Zosterops species across Kenya, including three species inhabiting highland sky-island systems; Z. kulalensis, Z. mbuluensis and Z. silvanus (one individual from stable habitat in Mt. Kasigau and one from a remnant forest fragment surrounded by disturbance in the Taita Hills, Fig. S2 location J at~1760m asl), and the Z. silvanus lowland congener, Z. avilateralis.
Genomic DNA was extracted from blood samples using a salt extraction protocol and whole genomes sequenced using paired-end Illumina on the NovaSeq 6000, generating high genome read coverage (>18x per individual). Raw reads were trimmed using cutadapt (v 2.10) and aligned to a Zosterops borbonicus reference genome (Leroy et al., 2019), using the BWA mem algorithm with default settings (Li & Durbin, 2009). Duplicate reads were marked using Picard MarkDuplicates (v 2.23.3). Mapping was checked using Qualimap (v 2.2.1). Resulting bam alignment les were used to generate consensus genome sequences for each individual (fastq) using the mpileup command in SAMtools and the vcf2fq command from vcfutils.pl, with the Zosterops borbonicus genome assembly as the reference. Each fastq le was ltered for sequencing errors by excluding sites at which the per-site root-mean-square mapping quality was below 25, the inferred consensus quality was below 20, and the variant read depth was less than 10x or more than twice the average across the genome. PSMC analyses were performed using the following xed parameters across each individual: maximum number of iterations (N) of 30, maximum coalescent time (t) of ve, initial theta/rho ratio (r) of one and parameter pattern (p) of`4+30 * 2+4+6+10'. These values were chosen in line with PSMC analyses conducted across other avian species (36 avian species; see Nadachowska-Brzyska et al., 2015). To scale the PSMC output to real time we use estimates of generation time as one year (as in main text) and neutral mutation rate as 0.2 × 10 −8 (Smeds et al., 2016) in the psmc_plot.pl command.

Results
Genomic inference of highland-lowland hybridization in isolated patches: The PSMC results provide some genomic evidence of demographic histories in highlandlowland White-eyes (Zosterops sp.). The lowland congener Z. avilateralis shows highest eective population size (N e ) with a peak at around 100 kya before decreasing. Independent of divergence times, sky-island specialists sharply lose N e and converge at comparably low levels (i.e. N e < 10 5 ). Interestingly, there is a separation of the otherwise highly congruent curve between Z. silvanus inhabiting the highly fragmented Taita Hills compared to the pristine Mt. Kasigau. The PSMC model suggests initial divergence of Z. silvanus inhabiting both the Taita Hills and Mt. Kasigau with low eective population size, while the lowland congener Z. avilateralis experienced rapid growth. In the more recent past, we observed an increase in eective population size in Z. silvanus in the Taita Hills, while all other highland species contracted, despite the most restricted geographical range being in the Taita Hills. This increase in N e could be possibly explained by a) range expansion or b) hybridization. While the former reason can be excluded given that the Taita Hills have been isolated for millions of years (i.e. there is no evidence that spread of the mountain cloud forest into the surrounding savannah habitat has occurred in the past), introgressive hybridization with another species might be the most plausible scenario. To this end, we expect introgressive hybridization, as a form of genetic rescue, from the lowland congener Z. avilateralis into Z. silvanus inhabiting these small forest fragments, an assumption that forms the logical basis for the scenario building reported in this study.